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Abstract 

The behavior of geodesic curves on even seemingly simple surfaces can be sur- 
prisingly complex. In this paper we use the Hamiltonian formulation of the 
geodesic equations to analyze their integrability properties. In particular, we 
examine the behavior of geodesies on surfaces defined by the spherical harmon- 
ics. Using the Morales-Ramis theorem and Kovacic algorithm we are able to 
prove that the geodesic equations on all surfaces defined by the sectoral har- 
monics are not integrable, and we use Poincare sections to demonstrate the 
breakdown of regular motion. 
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1. Introduction 

The study of the qualitative behavior of geodesies goes back to the begin- 
nings of differential geometry; in fact Clairault's relation dates from the 1740's. 
One of the many interesting features of the geodesic equations is their com- 
plexity: how can we solve such difficult nonlinear differential equations, and 
indeed what does it mean to 'solve' a differential equation at all? In the gen- 
eral development of Lagrangian/Hamiltonian dynamics the issue of integrability 
came to the fore, and the use of first integrals to solve a dynamical system via 
quadratures was developed. By integrability here we mean intebrability in the 
sense of Louiville: a Hamiltonian system with n degrees of freedom is Louiville 
integrable if there exist n independent first integrals of the motion in involu- 
tion (vanishing pairwise Poisson brackets). This formalism in turn lead to the 
consequences of non-integrability, the onset of chaos and the celebrated KAM 
theory. For a particularly interesting introduction to this topic see Berry 
McCauley [HI or Jose and Saletan (l7| . 

The integrability of the geodesic equations on some surfaces, and therefore 
the regularity of the geodesic trajectories, has been well understood for many 
years. Clairault's relation provides a second integral for surfaces of revolu- 
tion, and Jacobi [l8[ integrated the triaxial ellipsoid. These results have been 
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extended to rt-dimensional ellipsoids and to quadrics in general (see Audin Q), 
where interestingly degenerate cases can prove more difficult to analyse (see 



Davison et al [12[). Recently, the integrability of the geodesic equations has 



been demonstrated for more abstract manifolds, defined in terms of Lie groups 
(see for example Thimm |3l|). Aside from being interesting in their own right, 
the geodesies of a manifold also have some useful applications. For example, 
Maupertius's principle [7] relates the dynamics on an isoenergetic level of any 
mechanical system to the geodesies of a manifold, and the relativistic formula- 
tion of gravity dictates that massive particles move along the timelike geodesies 
of spacetime 



There are many powerful results regarding the non-integrability of the geodesic 
flow^ on manifolds of a certain class, for example a theorem of Kolokol'tsov (l9| 
says that two-dimensional compact surfaces of genus g > 1 do not admit geodesic 
flows that are integrable with a second integral polynomial with respect to the 
momenta (see Bolsinov [6] for a review and references). What's more the inte- 
grable geodesic flows on surfaces with g = 0, 1 with second integrals either linear 
or quadratic in the momenta are completely described. Nonetheless, new ex- 
amples of integrable geodesic flows with second integrals cubic in the momenta 
have recently been found (see for example Dullin and Matveev 13[) and the 
possibility of non-polynomial, and in particular meromorphic second integrals, 
is open. 

A very significant recent development with regards integrability has been 
the theorem of Morales Ruiz and Ramis 2|| 29[ . This theorem provides a link 
between the question of whether a nonlinear Hamiltonian system is integrable 
(with meromorphic first integrals) with whether the variational equations of a 
non-trivial solution to the Hamiltonian system are solvable. This is a great 
advance as the variational equations are linear with variable coefficients, and 
the solvability of equations of this type may be addressed using differential 
Galois theory. This approach has proved very successful (we give here only 
some examples) : in Morales Ruiz and Ramis [29[ the authors use this formalism 
to prove the non-integrability of certain problems from mechanics and celestial 
mechanics; in Morales Ruiz et al [3(| the non integrability of Hill's problem 
is shown; in Pujol et al [27[ the non-integrability of the swinging Atwood's 
machine is demonstrated; in Acosta-Humanez et al [1| some problems in celestial 
mechanics are studied and in Bardin et al 4] the generalized Jacobi problem 
is studied. In the majority of these papers, the solvability of the variational 
equations is determined using Kovacic's algorithm, and this is the approach 
we will follow in this work, with a description of the algorithm and relevant 
references given in Section 4. 

To the best of our knowledge this powerful Morales-Ramis formalism has 
not been applied to the study of geodesies on manifolds and their integrabil- 
ity properties, and the approach taken in this work is novel in this respect. 



1 Technically, the geodesic flow is defined on the cotangent bundle whose projection onto 
the manifold gives the geodesic curves, however we will use these two tems interchangeably. 
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In the following we will examine the properties of geodesies on surfaces de- 
fined in terms of the spherical harmonics, and we shall rigorously prove using 
the Morales-Ramis theorem the non-integrability of the geodesic equations on 
surfaces defined in terms of sectoral spherical harmonics. Thus we present a 
two-parameter family of surfaces which appear very simple and symmetrical 
and yet have complex and indeed chaotic geodesic curves on them; this is a 
novel addition as the surfaces for which non-integrability of the geodesic flow 
has been rigorously proved is small. Finally, we will also use Poincare sections 
to display the breakdown of regular motion as well as to understand the various 
families of closed geodesies on the surfaces; again this approach is novel in the 
case of geodesic flows. 

The layout of the paper shall be as follows: In Section 2 we review the 
formulation of the geodesic equations and describe the classes of surfaces which 
will be examined. In Section 3 we use Poincare sections to exhibit the breakdown 
in regular motion suggesting non-integrability, as well as describe the main 
families of closed geodesies on the surface. Finally in Section 4 we will review 
some necessary theory and prove the main result: the non-integrability of the 
geodesic equations on sectoral harmonic surfaces with n > 1. 

2. The geodesic equations 

One geometrical description of a geodesic curve is as follows: of all the curves 
connecting two nearby points on a manifold, the one that extremizes the length 
along that curve is the geodesic (minimizes in the case of a positive definite 
manifold). An alternative dynamical description is to imagine a test particle 
moving in the manifold free from external forces; if the particle's speed is con- 
stant, then the path it follows is a geodesic. Consequently, we may define a 
(dynamic) Lagrangian C in the following way: let the coordinates x a parame- 
terize the n-dimensional manifold (A4 , g) and let s parameterize a curve in the 
manifold, such that x a = x a (s). The length of the tangent vector to this curve 
is 



but this is simply the speed. Here g a b are the components of the metric tensor 
in the x a coordinate system, and from now on a dot denotes differentiation 
w.r.t. s and summation over repeated indices is implied. We may define the 
kinetic energy as one half the square of this quantity (assuming unit mass), and 
therefore the Lagrangian function is 



Putting the Lagrangian function into the Euler-Lagrange equations we arrive at 
the geodesic equations, 
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where the Christoffel symbols (metric connections) are given by 

r? c = \g ad {9bd,c + 9cd,b-9bc,d)- (3) 

The equivalent Hamiltonian formulation is clear, by simply letting % — C. For a 



more detailed derivation of the geodesic equations see Paternain 25], Do Carmo 
[Io| or alternatively Wald (33| . We note that some authors prefer to use the 
coordinates of the space the manifold is embedded in, and place a constraint on 
these coordinates (see for example Davison et al |l2|). This approach does not 
suit our purposes here. 

2.1. Surfaces defined in polar form 

Let us consider two-dimensional surfaces given in polar form, that is 

r = r(M), (4) 

where (9, </>) are the standard spherical coordinates with < 9 < tt, < 4> < 27r. 
We may calculate the geodesic equations in the following way: the line element 
for R 3 in spherical coordinates is given by 

ds 2 = dr 2 + r 2 d9 2 + r 2 sin 2 aAA ~ 2 



On the surface we have r = r(9, 4>) and hence dr = r_gd9 + r^d(f> and so 

ds 2 = (r 2 e + r 2 )d9 2 + 2r fi r^d9d<P + (r 2 + r 2 sin 2 9)d(j) 2 , (5) 

from which the metric coefficients may be read off, inserted into the Christoffel 
symbols in ([3]) which in turn are put into the geodesic equations in ^ , namely 

9 + T e J 2 + 2T e H n + T 9 ^ 2 = 0, (6) 

+ r^ 2 + 2r^0 + r^ 2 = o. (7) 

It will be most informative for later analysis to describe closed geodesies on 
these surfaces. For this, we note the following: if a smooth surface has a plane of 
symmetry, then the curve formed by the intersection of the plane of symmetry 
with the surface will be a geodesic. Geometrically the reason for this is clear: if 
the surface is smooth the plane of symmetry will be normal to the surface, and 
any normal curve is geodesic. For the sake of the analysis to follow however we 
will make the following argument: let us arrange the coordinate system so that 
the plane of symmetry is given by 9 — 7r/2; on this plane r # = 0. From (|6]), we 
see 9 — tt/2,9 — will be a solution if T^(6> = tt/2) = 0. Direct calculation 
from ((3} and (0 gives 



r w - 



r ,e( rr ,<P>P sin 2 9 — 2r 2 ^ sin 2 9 — r 2 sin 4 9) — cos 9(r 3 sin 3 9 + rr 2 ^ skU 
r(r 2 sin 2 9 + r 2 e sin 2 9 + r 2 , ) 



which vanishes on 9 = tt/2. This is equivalent to the 'invariant plane' con- 
struction typical in dynamics which we will describe more of below; we will 
subsequently refer to such geodesies as planar geodesies. 

Of course, to describe the geodesic completely we must solve the geodesic equa- 
tions reduced to this plane, namely (f> + ^i,(6 = 7r/2)0 2 — 0; we will address 
this issue in Section 4. 
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2.2. Surfaces defined by the spherical harmonics 

The surfaces we will examine in this work are those defined by the spherical 
harmonics in the following polar form 

r(M) = l + e*lm(M). (8) 

Here £ is a small parameter and Yi m is the spherical harmonic of degree I and 
order m (see Wang and Guo given by 



Ylm{e > 0) = J (2 4t(Z + m)7 )! P ^( cos0 y m ^ m = 0, ±1, . . . , ±/, 

where P™ are the associated Legendre polynomials. We may absorb the or- 
thonormality coefficient into e and take the real part of the exponential to write 

r((9,0) = 1 + eP[ n (cos 0)cos(m0), < £ < 1. (9) 

With a slight abuse of language, we shall refer to surfaces defined in terms of 
the spherical harmonic functions as 'spherical harmonic surfaces' and so on. We 
chose to study these surfaces for the following reasons: 

1. The spherical harmonics are complete, that is any function continuous over 
the sphere may be decomposed into a (possibly infinite) set of spherical 
harmonics. As such, general surfaces of a certain class can be decomposed 
into a number of spherical harmonic surfaces. 

2. The spherical harmonics naturally have the smoothness and regularity 
required for this analysis, and are non-singular as long as e < 1. 

3. The small parameter allows for a perturbative approach: with e = the 
surface is a sphere, and as we let £ grow we are deforming the surface away 
from the sphere. 

4. The spherical harmonic surfaces are compact, enabling the Poincare sec- 
tion analysis of the next section. 

There are three classes of spherical harmonic: the zonal, the sectoral and 
the tesseral (see Figure Q] for representative examples) . 

Zonal These are the spherical harmonics with m = 0, and thus r = 1 + 
e Ff(cos0). Since = 0, this means dC/dtfi = 0, which from the 
Euler- Lagrange equations leads to d£/d<j) =constant. This is of course 
Clairault's integral, as the zonal harmonic surfaces are surfaces of revolu- 
tion; the motion of the geodesies on zonal harmonic surfaces will therefore 
be regular. 

Sectoral These are the spherical harmonics for which I = m. We shall relabel 
the degree/order n in this case, and Y™ takes on a simple form: the sectoral 
harmonic surfaces are defined by 

r = 1 +e sin ,l (6»)cos(n0). (10) 
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Figure 1: Examples of spherical harmonic surfaces, for various values of e: on the left the 
I = 2, m = zonal harmonic, in the centre the I = m = 4 sectoral harmonic, and on the right 
the I = 7, m = 4 tesseral harmonic. 

This describes a two-parameter family of surfaces which admit planes of 
symmetry; as such they provide an interesting case for analysis and will 
form the main part of this work. 

Tesseral These are the remaining members of the family, with m = ±1, . . .,±1 — 
1. We will not examine these surfaces in any detail in this work, but discuss 
them briefly in the conclusions. 

The Hamiltonian function for the geodesic equations on the sectoral har- 
monic surfaces is given by 

m = geeO 2 + 2g e4> H + g^ 2 (11) 

where 

gee = (1 + ecos(?i0) sin™(6»)) 2 + e 2 n 2 cos 2 ^) sin 2 "~ 2 (6»), 

9e4> = g<t>e = — £ 2 n 2 cos(n^>) sin(r«/)) cos(6>) sin 2 ™" 1 (0), 

g^ = sin 2 (9) (1 + ecos(n0) sin"(6»)) 2 + e 2 n 2 sin 2 {n<f>) sin 2n {6). 

It is clear from (fTU|) , or indeed Figure [I] that the nth sectoral harmonic 
surface has n + 1 planes of symmetry: n vertical/meridianal planes given by 
<p = ni/n, i = 0, ...,n — 1 and 1 horizontal plane given by 9 = tt/2, the 
equator. Each of these planes defines a planar geodesic as described above. 

3. Poincare Sections 

This is a much used technique to attempt to visualize the whole of phase 
space in the case of a 2-degree of freedom Hamiltonian system. By fixing a 
value of % — Ho, & 3-dimensional submanifold of phasespace, and intersect- 
ing this manifold with a hyperplane (typically given by a coordinate taking 
on a fixed value), we consider successive intersections of solution trajectories 
with this hyperplane. If the Hamiltonian system is integrable, then the second 
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integral restricts the solutions further to a 2-dimensional submanifold of phas- 
espace, whose intersections with the hyperplane will trace out a closed curve. 
On the other hand, if the Hamiltonian system is not integrable then successive 
intersections of the solution trajectory with the hyperplane may instead fill a 
2-dimensional region. 

Typically, for a fixed set of parameters a number of Poincare sections for 
various values of the 'energy' T-Lo must be drawn, for example in the Henon- 
Heiles system [l6j]. This is not necessary in the current situation, since the 
numerical value of the Hamiltonian function may be changed from "Ho to "Hi by 
the change of independent variable s — > ^/'Ho/'Hi s. Therefore the numerical 
value of H signifies the parameterization of the geodesic curve, rather than any 
'energy' concept. Indeed, any affine coordinate transformation s — > as + (3 will 
change the value of H but leaves ([2]) invariant. Therefore, any qualitative change 
in the Poincare sections will be due to variations in the surface's parameters {n 
and e for the sectoral harmonic surfaces), rather than the value of H. In what 
follows, we will assume the geodesic is parameterized by arc-length (sometimes 
referred to as 'unit-speed curves' in the literature), and therefore 

U=\. 

In choosing the hyperplane to act as the Poincare section we must ensure 
all solutions intersect the hyperplane trans versally; no tangential grazings are 
allowed, and a guarantee that the trajectory will intersect the surface again is 
preferable. Any of the n + 1 planes of symmetry would be suitable candidates, 
although for simplicity we will use the equatorial plane (some comments regard- 
ing the use of mcridianal planes are given in Appendix A). First we prove the 
following lemma: 

Lemma 1. Along a geodesic, 9 cannot obtain a maximum/minimum in the 
northern/southern hemisphere respectively, at least for moderate values ofe. 

Proof. We remind the reader that 9 is measured down from the positive z- 
axis, and the northern and southern hemispheres are given by 9 € (0,7r/2) or 
(7r/2, 7r) respectively. Due to symmetry, we need only show that 9 cannot obtain 
a maximum for 9 6 (0,7r/2). 

From @, a maximum w.r.t. 9 in the northern hemisphere would require F^ > 
for some 9 <E (0, n/2) (on the sphere, Ti i — - cos 6* sin 6» < for < 9 < n/2, 
ruling out this possibility). We can write 

g rcos0sin 3 6) . 

= del(g) A sm ^ C0S W > ); n ' £ )- 

As the surface is positive definite to prove the lemma we must show / > in 
O = {0 < 9 < 7r/2, < 4> < 27r}. Letting p = sm9 with < p < 1 in O, we can 
write 

f(p) = 1 + a lP n + a 2 p 2n - 2 + a 3 p 2n + a^ 2 + a 5 p 3n 

where the depend on cos(n</>), n and s ((/) here being the coordinate at which 
the proposed maximum takes place). When cos(n0) > 0, all the ai are positive 
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and so f(p) > in f2. 

If, on the other hand, 4> is such that cos(n<^) < 0, Descarte's rule of signs tells 
us that f(p) can only have 1 or 3 roots in < p < 1, and since /(0) > this 
is impossible unless /(l) < 0. Letting /(l) = gives a cubic in cos(n0), and 
we can then find explicitly the value of e as a function of n at which the value 
of /(l) becomes negative for some values of (j>. If we call this critical value £*, 
then for example when n = 2, 3, 4 we find e* = 0.57, 0.497, 0.445. □ 

Note this lemma rules out the possibility of closed geodesies which do not in- 
tersect the equator at some stage (the only 9 = constant geodesic is the equator) . 
What it does not rule out, however, is the possibility of a geodesic asymptot- 
ically approaching the equator. In fact, this cannot occur as the equator is of 
elliptic type. A rigorous proof of the non-hyperbolic nature of the equator is too 
much of a digression; certainly numerical evidence can be easily provided using 
Floquet theory (see [35J for a detailed description) by calculating the eigenval- 
ues of the monodromy matrix and showing they remain on the unit circle in the 
complex plane for moderate values of e > and various values of n. Having said 
that the most immediate demonstration of the elliptic nature of the equator is 
to be found in the Poincare section shown in Fig. A. 3. 

Proposition 1. With the equatorial geodesic of elliptic type, any geodesic which 
intersects the equatorial plane transversally must re-intersect the equatorial plane 
an infinite number of times, at least for moderate values of e. 

Proof. First we can see that any trajectory which intersects the equatorial plane 
transversally cannot subsequently graze this plane tangentially as it is an invari- 
ant plane. If we consider a geodesic which crosses the equator with 9 < (i.e. en- 
tering the northern hemisphere) then there are two possibilities: 1) the geodesic 
does not pass through the north pole, and 2) the geodesic passes through the 
north pole. In the event of case 1), the geodesic must attain a minimum w.r.t. 
9 {9 is an initially positive decreasing continuous function of arc-length and so 
must either: attain a minimum; pass through zero (case 2); or asymptote to 
some constant value which is impossible as the surface is compact and there 
are no equilibrium points or closed geodesies contained in a hemisphere). As 
the geodesic cannot attain a maximum w.r.t. 9 in the northern hemisphere or 
asymptote to the equator, it must subsequently re-intersect the equator. In 
the event of case 2), consider a geodesic which starts at the north pole. While 
there is a degeneracy in defining 4> at the poles there is no problem in defining 
4>. Following the geodesic in this direction it must intersect the equator for the 
reasons given in case 1), but the same can be said if we follow the geodesic 
"backwards" , and so it must intersect the equator in both directions. □ 

A final consideration is the ranges of velocities the geodesies may explore. 
On any of the planes of symmetry, either r t $ or r.^ vanishes, thus (from ([5|)) 
getj, = 0. Therefore on the section we have 

gee9 2 +W> 2 = 1, (12) 
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which is an ellipse of semi- axes l/^/g$g and \j Since g^(B — tt/2) = 
(1+e cos(n0)) 2 +e 2 n 2 sin 2 (n^>), when choosing our initial conditions to generate 
the Poincare section plots we should therefore choose random values for <j> in 
the range 4> G {~4>m, 4>m) with 

9m = < i 

/ „2 (1 + e 2 ( „2_7 

The initial conditions can now be chosen in the following way: random values 
for are chosen in the range (0, 2ir), random values for <p in the range (— </> m , m ), 
and the value of 9 is found from % = \ (we will always take the positive value to 
give a 'downward' pointing geodesic). The geodesic with these initial conditions 
is followed for a long time and intersections with the equatorial plane with 6 > 
are monitored. We show in Figure [2] some of the resulting Poincare sections for 
the case of the n — 3 sectoral harmonic surface. Some interesting points to note 
are: 

1. There are three main families of closed geodesies; these appear as fixed 
points of the Poincare section, and are as follows: Firstly there are n 
planar geodesies as described in Section 2 leading to fixed points at <j) = 
ni/n, i = 1, . . . ,n and <f> — 0. Secondly on the sectoral surfaces with n odd 
there are n perpendicular geodesies, which are not planar but intersect the 
equator at right angles, leading to fixed points roughly midway between 
those due to the planar geodesies. Thirdly there are oblique geodesies 
which are simple (period- 1) and give rise to fixed points off the horizontal 
axis of the Poincare section. 

Of course, there are an infinite number of closed geodesies on topological 
spheres (see Bangert [3j and Franks [15j ) and we make no attempt to 
classify them all. 

2. The stability of these closed geodesies is perhaps not what would be ex- 
pected (Note: we are using only the Poincare sections to draw conclusions 
about stability; a more rigorous stability analysis would make an interest- 
ing topic for future research). Consider first the sectoral surfaces with n 
even, as shown in Figure 1 (b). The planar geodesies alternate between 
those for which d 2 r/d<fi 2 < and those for which d 2 r/dtfi 2 > every- 
where along the curve. In analogy with the surfaces of revolution [lOj], 
we might therefore expect that these geodesies would alternate between 
stability and instability, respectively. However this is not the case: all 
planar geodesies are in fact stable, leading to elliptic fixed points on the 
Poincare section. 

If we now consider the sectoral surfaces for n odd, then d 2 r/dcj> 2 is both 
positive and negative along each planar geodesic. In fact, these planar 
geodesies are all unstable; what's more each is homoclinic to itself (at 
least for small e, see below) and thus admits two elliptic points within its 
homoclinic loops, see Figure [2] These elliptic points form a fourth family 



if s < 



if e > -4- 



(13) 



n 2 -f 
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of closed geodesies, and some numerical experimentation reveals these are 
also period- 1. 

The perpendicular geodesies mentioned in 1. are all stable, leading to ellip- 
tic fixed points in the Poincare section, and most importantly the oblique 
geodesies are unstable, leading to hyperbolic fixed points. These hyper- 
bolic fixed points are heteroclinic to one another as Figure [5] shows (in 
fact there are complex heteroclinic networks, where up to 2n saddles are 
connected to one another), and this is the crucial ingredient leading to 
chaotic motion, which we shall describe next. 

As is clear from Figure^ chaotic motion takes place for large e, that is 
as we deform the surface away from the sphere. The onset of chaos as 
depicted in Figure [5] is straight out of a textbook: the invariant manifolds 
of unstable fixed points intersect leading to a heteroclinic tangle, and as e 
increases further rational invariant tori are destroyed. This phenomenon 



is well described in the literature, see for example Jose and Saletan 17] . 
Calkin or Berry @. The value of e at which this onset of chaotic motion 
takes place varies for different values of n, but note that we see chaotic 
motion for values of e less than those mentioned in Lemma 1. 



The onset of chaotic motion is visible in the Poincare section of all sectoral 
harmonic surfaces, with the exception of the n = 1 surface. On this surface the 
geodesic flow is perfectly regular for all values of e € [0,1). This is because the 
n = 1 sectoral surface is in fact a surface of revolution; it is the surface formed 
by revolving the limacon around the cc-axis, as we may show: 

The n = 1 sectoral surface develops a 'dimple' along the x-axis. Let us rotate 
the axes (or alternatively the surface) about the y-axis so that the dimple is along 
the 2-axis. As this rotation amounts to (x,y,z) — > (z,y,—x) we can define a 
new set of spherical coordinates (9, cj)) related to the old spherical coordinates 

by 



sin 6 = y 1 — sin 2 9 cos 2 0, cos <fi = _ - . 

V 1 — sin 2 9 cos 2 <j> 

The n — 1 sectoral harmonic surface is therefore given by 

r = 1 + e sin 9 cos <f> = 1 + e cos 9, (14) 

the parametric form of a limacpn, and clearly dr/dtp = 0. This surface will 
therefore have integrable geodesic flow, a fact that will prove useful in the next 
section as a check of the procedure carried out. 

The Poincare sections provide strong evidence for the following conclusion: 
the geodesic equations for the sectoral harmonic surfaces with n > 1 are not 
integrable. However, the Poincare sections cannot be used as proof of this 
statement, for the following reasons: 

1) as the break up of regular motion apparently happens only for larger values 
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(c) 

Figure 2: Sample equatorial Poincare sections for the n = 3 sectoral harmonic surface, with 
the following values of e: (a) e = 0.1; (b) e = 0.2; (c) e = 0.3. The axes are <f> and <f>. See text 
for discussion. 
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of e, there are possibly small values of e for which the equations are integrable, 
as in the general Henon-Heiles system for example [23 |; 

2) the integrations which generate the plots are run over a hnite time; if the 
integrations were run over a longer time interval perhaps a sequence of points 
which seems to form a closed curve would begin to wander away. This 'stickiness' 
phenomenon is in fact quite common in chaotic systems (see for example Murray 
and Dermott [3]), and 

3) while great pain may be taken to ensure accuracy in numerical investigations, 
there is no escaping the fact that these are only approximate solutions to a 
system of nonlinear ODE's. 

Taking these points on board we seek a rigorous and analytical test of the 
integrability or otherwise of the geodesic equations. This will be the goal of the 
next section. 



4. The variational equations and non-integrability 

While proving a dynamical system is integrable is very challenging, and 
amounts to essentially finding the first integrals, or showing the Hamiltonian is 



separable in some sense (see for example O Mathuna 2l|), there are a number 
of ways of proving a dynamical system is not integrable. This is because in- 
tegrable systems are very structured, and must admit a number of properties; 
showing one of these properties does not hold suffices to prove the system is not 
integrable. While Poincare famously showed the non-integrability of the 3-body 
problem over 100 years ago [2a |. a particularly fruitful modern approach is that 
due to Morales and Ramis, as exp ounded in (23[,[2i| and (29[, and applied to 



the 3-body problem in [8[ and [32j among others. 

This method focuses on the variational equations which arise from linearizing 
around a non-constant solution to the dynamical system. In loose terms (which 
we will define more precisely below), if a Hamiltonian system is integrable then 
the variational equations will be solvable, in the sense of differential Galois 
theory. Consequently, if the variational equations are not solvable then the 
Hamiltonian system is not integrable. This means we need only focus on the 
variational equations, which is a great advantage as the variational equations 
are of course linear and therefore much easier to analyse. We will give a brief 
summary of the required results first. 

4-1. Solvability and the Morales-Ramis theorem 

Consider the system of linear ordinary differential equations 

f = A£, ' = d/dz (15) 

where the elements of A are in some field K (typically the field of meromorphic 
functions, C(z)). We may construct the Picard-Vessiot extension to this field 
L/K by adjoining to K the solutions to (|15p. Associated with (|15|l is the 
differential Galois group Q which is the group of automorphisms of L / K which 
leave fixed the elements of K . The linear system is solvable if we may go from 
the base field K to the field extension L/K by adjoining to K 
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1. integrals, 

2. exponentiation of integrals, and 

3. algebraic functions 

of the elements of K ([23[ gives a complete introduction to this topic). Thus we 
can 'build' the solutions to (TT51) from the elements of the field that the coeffi- 
cients of (|15p belong to, by performing a finite number of 'allowed' operations. 
Solutions of this type may be referred to as closed form or Louivillian solutions, 
and the differential equation may be called solvable or Louivillian. 

We define the variational equations as follows: for the dynamical system 
dx/dt = f(x) with non-constant solution 7, we let x = 7 + £ and the first 
variational equations are given by 

7 



dt V dx 



The variational equations can be separated into a tangential and a normal com- 
ponent; as autonomous Hamiltonian systems admit the Hamiltonian itself as a 
first integral this is 'inherited' by the tangential variational equation which is 
therefore always solvable. And now importantly, the normal variational equa- 
tion will inherit any other first integrals, if they exist, and will therefore also be 
solvable. This is the Morales- Ramis theorem, which we will quote from [23]: 

Theorem (Morales- Ramis). For a In dimensional Hamiltonian system as- 
sume there are n first integrals which are meromorphic, in involution and inde- 
pendent in the neighborhood of some non-constant solution. Then the identity 
component of the differential Galois group of the normal variational equation is 
an abelian subgroup of the symplectic group. 

In other words, the normal variational equation (NVE) is solvable in the 
sense described above. One advantage of the MR theorem is that it considers 
the differential Galois group Q associated with the differential equation (|15p. 
which is a Lie group, rather than the monodromy group used in Ziglin's theorem, 
which is discontinuous. 

This theorem suggests a procedure for testing the integrability of a Hamil- 
tonian system: find a non-constant solution; linearize around it and separate 
out the NVE; test the NVE for solvability. If the NVE is not solvable, then the 
Hamiltonian system is not integrable. To test the solvability of the NVE, we 
will use the Kovacic algorithm. 

4-2. Kovacic algorithm 

The Kovacic algorithm can be used to find the closed form solutions of a 
second order linear ODE with rational coefficients, and is complete in the sense 
that if the algorithm does not find such solutions then none exist. The original 
source for the algorithm is Kovacic [2(j , but see also 14 [ , [2^] , Q , [3 , to name 



a few, for presentation, discussion, and application of the algorithm. As we will 
see, the NVE associated with the system studied in this work will be Fuchsian, 
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which allows for a more concise version of the algorithm. As such, we will follow 
most closely the presentation of the Kovacic algorithm as given in Churchill and 
Rod [ll|, and note that the algorithm is also valid in the non-Fuchsian case. 
Consider the linear differential equation 

S" = r(z)S, reC(4 (17) 
If this equation is Fuchsian we can write 



^ («-Oj) 2 



where k is the number of finite regular singular points at locations z = ctj. When 
^2 Sj = then z = oo is also a regular singular point, with = + $j a j)- 

The indicial exponents are 



at z = <Zj and z = oo respectively. Kovacic [20| proved the following theorem: 

Theorem (Kovacic). Let Q be the differential Galois group associated with 
(|17[) . and note Q € SX(2, C). TTien on/?/ one of four cases can hold: 

1. is triangulisable (or reducible), in which case (fTTJ) /ias a solution of the 
form ^ — e^ with lu £ C(z) (we could say u> is algebraic over C(z) of 
degree 1 ); 

2. Q is conjugate to a subgroup of 



A" 1 /' ' 1 J j 

in which case (|17p /ias a solution of the form £ = where lu is algebraic 

over C(z) o/ degree 2; 
3. is finite, in which case (1171) /ias a solution of the form £ = e-f w where 

U) is algebraic over C(z) o/ degree 4, 6 or 12; 
A. Q = SL(2,C), in which case (|17| is not solvable. 

To show the NVE is not solvable, we must check that each of the cases 1-3 
cannot hold. If this is the case, then Q = SX(2,C) whose identity component 
(also SX(2,C)) is not abelian, and therefore the original Hamiltonian system is 
not meromorphically integrable according to the Moralcs-Ramis theorem. 

Kovacic provides an algorithm whereby each case is checked in two stages: 
first, we consider certain combinations of the indicial exponents and retain all 
those leading to a non- negative integer d, i.e. d £ No- Second, for each of these 
values of d we attempt to construct a monic polynomial of degree d, or possibly 
a set of polynomials, which satisfy certain equations. If successful, the solution 
to (fT7|) can be found. We will state the details for each case in turn: 
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Case 1 Define the modified indicial exponents 

af = rf if ft ^ 0; af = 1 if ft = and 8j ^ 0; af = if ft = ft = 0, 
°4 = T ta ^ A» 7^ 0; a+ = 1, a~ = if ft^ = 0. 
Find all combinations of these modified indicial exponents such that 

k 

d = a to - a t 

is a non-negative integer. For each d found, look for a monic degree d polynomial 
P(z) satisfying 



P" + 26P' + (9' + 9 2 - r)P = 0, 9 = 




If found, a solution to ((171) will be £ = e ^ " where uj = 6 + P'/P. 
Case 2 Define the following sets: 

E 3 = {2 + ey/1 + 4ft, e = 0, ±2} n Z if ft 7^ 0; 

E 3 = {4} if ft = 0, ft 7^ 0; Ej = {0} if ft = ft = 0; 

Eoo = {2 + eVl + ^oo, e = 0, ±2} n Z if Poo + 0; 
Soo -{0,2,4} if Poo = 0. 

Find all combinations of ej G -Ej and e,^ € i?oo, not all even integers, so that 

1 k 
d= 7;{ e °°~J2 e j) 
3=1 

is a non-negative integer. For each d found, look for a monic degree d polynomial 
P(z) satisfying 

1 k 

P"'+WP"+(W 2 +?>6'-Ar)P'+(0"+ , Z66'+6 :i -Ar0-2r')P = 0, 6 = - V — ^— 

2 ^— ' 2 — a, 

3=1 J 

If found, a solution to p7|) will be £ = e-f w where a; is a root of uj 2 + 4>uj + (^<fi' + 
\(j> 2 -r) = Oand<£ = + P'/P. 

Case 3 In this case we must follow the procedure below for each of the subcases 
N = 4, 6, 12. Define the following sets: 

Fj ={6 + ^Vr+4ft,e = 0,±l,...,±f}nZif ft 7^0; 

Fj = {12} if ft = 0, ft ^ 0; Fj = {0} if ft = ft = 0; (19) 

Foo = {6 + ^Vl + 4ftx>,e = 0,±l,...,±f }nZ. 
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Find all combinations of fj € Fj and /oo G P^ so that 

3=1 

is a non-negative integer. Now define 

j=i 3 i=i 

Let Pat = — P where P(z) is a monic polynomial of degree d, and recursively 
define a set of polynomials P/v-i, Pzv— 2) • • • , P—i in the following way: 

Pi-! = -SPl + [(N - i)S' - S9]Pi - (N - + l)P+i, * = N, N - 1, . . . ,0. 

If we find P_i = 0, then use the set of polynomials Pi as coefficients in the 
following polynomial and find a solution w. 

A solution to ^T7J) will be ^ = e/ w . 

If these three cases fail to produce a closed form solution, it is because such 
solutions do not exist. We are now in a position to derive the NVE and apply 
the Kovacic algorithm to it. 

4-3. The normal variational equation 

First we must derive a non-constant solution to linearize around (there are no 
constant solutions to geodesic equations). In general, it is a difficult task to find 
an exact solution to geodesic equations, even for simple symmetrical manifolds. 
In the case of the sectoral harmonic surfaces however, we are greatly aided by 
the existence of certain planes of symmetry, leading to planar geodesies. These 
planes of symmetry play the same role in this context as invariant planes play 
in more general contexts. In what follows we will take the equatorial planar 
geodesic as our nominal solution (see Appendix A for some comments on the 
use of meridianal planar geodesies). 

In the equatorial plane 9 = ^,9 = 0, the geodesic equation is (from ©,([7]) 
and noting that = 0, where we use a tilde to denote the value on 9 = f ) 

<^ + f^ 2 = (20) 

which admits Hamiltonian 

2H = g^ 2 = 1. (21) 
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To solve for this equatorial planar geodesic we would need to rewrite (|2"Tj) as 



J ds = J d<pyfg^ = J d(j)\J (1 + scos{ncf))) 2 + e 2 n 2 sin 2 (n<?(>) (22) 

and solve to find s = s (</>), then invert this relationship to find cf> = <fi(s). 
The integrals in (|22|) can be written in terms of elliptic integrals of various 
kinds, at least for some values of n, however the analysis would benefit from 
a unified treatment with n as a parameter. That the integrals in (|22l) are 
daunting should come as no surprise: we are essentially trying to find the arc- 
length parameterization of the curve r = 1 + e cos(n0), which we could describe 
as a limacon with n 'lobes', and there are few curves for which the arc-length 
is known exactly. 

Fortunately we do not need to solve (j2"2")l exactly; for the sake of the NVE 
we can make a change of independent variable which requires then only (|21[) , as 
we will show. Let us denote the planar geodesic by 

(M,M) = (i,^o,0,<£o), (23) 

and let = 4> = 4>o + r). Linearizing around this solution we find the equa- 

tion in £ naturally decouples from that of i), and since the equatorial/invariant 
plane is given by 9 = J the variation normal to this plane is given by £, and 
thus the NVE is 

'(= (-P^ 2 ) £+(-2fg^) t (24) 

where the coefficients are evaluated along the nominal equatorial geodesic. Now 
we make the change of independent variable z — £cos(n</>(s)) so that 

d . . .d<f> d ensminif) d , . 

— = — ensm(ncp) — — = — — — (25) 

as ds dz ^/arhrh dz 



and so on. We may replace instances of cos(n0) with z/e, to give an NVE 
with rational polynomial coefficients (this process is sometimes referred to as 
algebrization, and preserves the identity component of the differential Galois 
group [3(|), which we may write as (letting £' = d^/dz etc.) 

e+p(zK' + q(^ = Q- (26) 

Finally, we transform this equation into the standard form required for the 
Kovacic algorithm 

£" = r(z)t 



where 



-q + jP 2 + y 



and we extend z to the complex domain. We are now in a position to state the 
main result of his paper: 
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Theorem 1. The geodesic equations for sectoral harmonic surfaces with n > 1 
are not integrable. 



Proof. When n > 1, the NVE has 6 regular singular points: five finite regular 
singular points and one at infinity. The finite regular singular points are 

{oi, ... ,a 5 } = {-1, e,-e,p+,p-} (27) 

with 

i± n yi + £ y-i) 

p+ '- n^—i ■ (28) 

The case n = 1 is special, and shall be dealt with below. When n > 1 and 
< e < 1 all singular points are real and distinct. The coefficients in the partial 
fraction expansion (|18|) are 



L -3 -3 5 5 I „ n+1 
ft = {°>T6'l6'16'16}' P-'-nT' (29) 

and the S^s are given in Appendix B; here we need only mention that ^ Si = 0. 

We must now check all three cases of the Kovacic algorithm to see is the 
NVE solvable. Often in the literature some of the cases can be ruled out straight 
away, based on the "necessary conditions" given in Section 2.1 of Kovacic [201 ] - 
However all of these necessary conditions are satisfied in the present case since 
the order of each pole of r(z) is 2 and moreover 

v^+4ft = {i,i \,\,\\ eQ, 71+^ = — eQ; (30) 

villi) n 

therefore we must examine each case. 

The first stage in each case consists of finding all combinations of the /3 
coefficients such that d is a non-negative integer. Unfortunately, there are a lot 
of such combinations. We show the combinations in Table 1, where N = 1 and 
N = 2 represent cases 1 and 2 respectively, and N = 4, 6, 12 represent case 3. 
The values of n, the specific sectoral harmonic under consideration, are shown 
in the column on the left. Each entry in the table gives the possible values of 
d in bold, and in brackets the number of combinations of indicial exponents 
which lead to this value of d. As can be seen immediately from the table, when 
n = 7, 8, 9, 11 and n > 12 there are no combinations leading to d £ No; as such 
the NVE in these cases is not solvable. 

For the other values of n we must try to construct a monic polynomial of 
degree d which solves the required equations. Clearly there are too many to go 
through individually; for the sake of demonstration we will describe the most 
involved case, n — 2, N — 12. The relevant sets defined in (fll?)) are 

F 1 = {12}, ^=^ = {6,7,5,8,4,9,3}, ^ = ^={6,9,3,12,0,15,-3}, 

Foo = {6, 8, 4, 10, 2, 12, 0, 14, -2, 16, -4, 18, -6} (31) 
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It is easy to see how there can be many combinations leading to d € No- The 
largest value comes from the combination 

d = 18- (12 + 3 + 3-3 - 3) = 6. 

In this case, we let P = z 6 + p 5 z 5 + . . . + p , and define 

P12 = -P, 

Pn = —SP[ 2 — S9Pi2, 
P 10 = -SPi, + [S' - S0]P n - \2S 2 rP 12 , 
P 9 = ... 

all the way down to P_i. These computations quickly become very large and 
cumbersome. A useful note is that each Pi is simply a polynomial in z, albeit 
one with very complicated coefficients. The order increases in a regular fashion: 
since S = 0(z k ), S9 — 0(z k ^ 1 ) and S 2 r = 0{z 2<<k ~ v> ) where k is the number 
of finite regular singular points, we see that Pn = 0(z d ), P/v-i = 0(z k ~ 1+d ) 
and so on until P_i = 0(z^ N+1 ^ k ~ 1 ^ +d ). In the present case this means P_i = 
0(z 58 ), although the leading coefficient always vanishes. Another useful point 
is that we only need the first few coefficients in P_i; setting the first equal to 
zero gives a value for p$, the second gives a value for p4 and so on until we reach 
a term which cannot be set equal to zero. Then we know that P_i ^= 0, which 
means that this combination for d does not produce a solution. 

This procedure has been followed for all combinations listed in Table 1 (us- 
ing Mathematica). For each and every case, it is not possible to construct the 
required polynomial P. As such, the NVE arising from linearising about an 
equatorial geodesic on the sectoral harmonic surfaces with n > 1 does not have 
a differential Galois group belonging to cases 1,2 or 3 in Kovacic's theorem; 
instead, Q = SX(2,C) whose identity component is not abelian. Therefore ac- 
cording to the Moralcs-Ramis theorem the geodesic equations are not integrable. 
This ends the proof. □ 

We can now deal with the special case n = 1: 

Theorem 2. The equatorial NVE on the n = 1 sectoral harmonic surface is 
solvable. 

Proof. When n = 1 the equatorial NVE has five regular singular points, four 
finite at aj = {— l,e, — e, p} with p = — ^(1 + e 2 ), and one at infinity. We 
calculate 

and hence 

^ = {(l,l),(f,I),(f,I),(f,-I),}, (33) 
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n 


N = 1 


N = 2 


N = 4 


N = 6 


N = 12 


2 


0(4) 


0(3), 1(1) 


0(4), 1(2), 
2(1) 


0(21), 1(10), 
2(3), 3(1) 


0(31), 1(20), 2(13), 
3(8), 4(4), 5(2), 6(1) 


3 








0(9), 1(3), 2(1) 


0(9), 1(6), 2(3), 
3(2), 4(1) 


4 




0(2) 


0(2), 1(1) 


0(7), 1(2) 


0(7), 1(3), 2(2), 
3(1) 


5 










0(2), 1(1) 


6 








0(3), 1(1) 


0(3), 1(2), 2(1) 


7 

8 
9 












10 










0(1) 


11 












12 








0(2) 


0(2), 1(1) 



Table 1: The possible values of the non- negative integer d; see text for an explanation. 



There is a combination giving 

d =!-( i+ j + i-i)=» < M > 

and since 8 thus defined solves the Riccatti equation 8' + 8 2 — r = we have 
the pair of Louivillian solutions 

Ci = C 1 (z + l)(z 2 -e 2 ) 3 / 4 (z- P )- 1 / 4 , ^ = C 2 ^J^ 2 dz. (35) 

This is of course entirely as expected since the n = 1 sectoral surface is a surface 
of revolution and thus has integrable geodesic flow. 

□ 

5. Conclusions 

The application of Moralcs-Ramis theory and Poincare sections, two tech- 
niques most often seen in the realm of celestial mechanics and mechanical sys- 
tems, to study the integrability of geodesic flow presents many interesting pos- 
sibilities. By choosing a specific class of surfaces to analyze we can rigorously 
prove the non-integrability of the geodesic equations. It would be of interest 
to extend the results of this work to examine further surfaces; for example the 
tesseral harmonics mentioned in Section 2 again have (meridianal) planes of 
symmetry, and the NVE will have coefficients that contain Legendre polynomi- 
als, and in general, it would be very informative if we could examine arbitrary 
surfaces with a plane of symmetry and try and arrive at some obstructions to 
integrability. We note that the presence of a plane of symmetry makes the 
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decoupling of the NVE straightforward, however Morales Ruiz [23] outlines a 
procedure in the absence of an invariant plane. The completeness of the spheri- 
cal harmonics means that surfaces of a very broad class can be decomposed into 
a set of spherical harmonics, and perhaps by analyzing these individually we 
can make some conclusions about their sum. Care must be taken however: the 
triaxial ellipsoid can be decomposed into a set of (infinite) spherical harmon- 
ics, each of which would likely have non-integrable geodesic flow, and yet their 
recombination gives a surface which is integrable. A further issue to consider 
is how the surface in question is described; the parametric description of the 
harmonic surfaces led easily to the calculation of the geodesic equations and in 
particular the separation of the normal variational equation. If the surface were 
defined in algebraic form, as is more suited to the double torus for example, this 
separation may be more problematic. Finally, this subject provides an inter- 
esting and fruitful combination of techniques from dynamical systems theory, 
differential geometry and differential Galois theory. 
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Appendix A. Meridianal sections 

The analysis above would carry through were we to use one of the meridianal 
planes instead of the equatorial plane, with two points to note: 

1. When drawing the Poincare sections, the fact that (j> =constant is only a 
half plane complicates matters. Better to rotate the coordinate system (or 
surface) to make the meridian an equator. For this we need to rotate the 
surface about the a;-axis, and for this we define a new set of coordinates 

&,<P by 

cos $ sin $ cos ip 

sin 6 



srnO = y 1 — sin $sin <p, sin0=— : — — , cos^> = 



Now to write the nth sectoral harmonic we expand cos(n0) and then make 
the replacement; for example, the 3rd sectoral harmonic changes as 

sin 3 6 cos(30) = sin 3 #(cos 3 0—3 cos (fism 2 (j>) = sin 3 $cos 3 ip — 3 sin ■& cos <p cos 2 

An example of a Poincare section for this rotated surface is shown below; 
the same families of closed geodesies can be identified. 
2. In deriving the variational equations we may linearize around a meridianal 
planar geodesic, to derive a normal variational equation of the form given 
in (|26[) . In the case of an equatorial planar geodesic the denominators 
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Figure A. 3: Sample Poincare section for the rotated n = 3 sectoral harmonic surface with 
e = 0.1. The elliptic nature of the equatorial geodesic (the fixed point at the centre of the 
plot) can be clearly seen. 



of the coefficients are at most quadratic in z, however in the case of a 
meridianal planar geodesic the denominators are polynomials of degree 
2n + 2. While the singular points are still regular, the number, location 
and behaviour of the set of singular points will vary with n, complicating 
matters. 



Appendix B. Coefficients of equation (|18|) 



1 ~ n(£ 2 -l)' 

_ 8 + n 2 ± 2e(8 + 4n + in 2 ) + e 2 (8 + 8n + 5n 2 + 8ri 3 + 4n 4 ) 
2/3 ~ 16e(l±e) 2 n 2 ' 

h ' 5 = - 2n(l-e 2 ) 2 ( £2(1 ~ 2 " 2) ~ 1 ± + c 2 Wl + e 2 (n 2 -l)] 
4 111 



«4/5 _ a i «4/5 — 0-2 a 4 / 5 — a 3 a 4 / 5 — a 5 / 4 



where the a,j are as defined in (|27]l .([28 ]l . 
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